N85 -2249 2 

SELF-CONSISTENT SIMULATION QC PLASMA INTERACTIONS WITH - 
SECQNDARY-EMITTlNl INSULATORS 


S-* T. Brandon, R. L. Rusk, T.P. Armstrong, and J. Enoch 
University of Kansas 
Lawrence, Kansas 66045 


A cylindrical particle-in-cell (PIC) plasma simulation code applicable 
to plasma- densities encountered in low earth orbit (LEO) is nearing 
completion. The simulated geometries include that of a plain disk and a 
disk surrounded by a. dielectric. Both configurations are mounted upon a 
ground- plate in contact with a plasma environment. Techniques allowing 
simulation of dielectric charging using PIC time scales are discussed. 
Current versus voltage characteristic curves are calculated and the results 
are compared to experimental data. 


INTRODUCTION 


The plasma densities present in LEO (10 - 10 /cn5 may cause the 
collection of large plasma coupLing currents for spacecraft operating at high 
voltages. In particular, large high-voltage solar- cell arrays exposed to the 
LEO environment could collect enough parasitic current from the ambient 
plasma to degrade their performance. Additionally, exposed, dielectric 
material can develop large differential potentials. Punctures in insulators 
over conductors have been seen to collect currents much larger than would be 
expected based on the area of the exposed conductor (ref. 1). 

In this paper, a status report of the results of a continuing effort to 
develop a self-consistent numerical simulation to explore more thoroughly the 
interactions of an ambient plasma with a conducting disk, which may be 
partially covered with a dielectric material, is presented. The disk 
surrounded by dielectric material represents a hole in an insulator covering 
a conductor. The background plasma parameters are chosen to resemble 
conditions of LEO. Since all particle trajectories are known, any process 
which can be modeled statistically for a single plasma particle can be 
included by the simulation code. Plasma interactions currently considered 
include the effects of charge collection on the dielectric and secondary 
electron emission. The simulation region includes a ground plate on which 
the disk or disk and dielectric are mounted. The resulting configuration 
closely matches the experiments (ref. 2). 

This combination of geometry and plasma parameters incorporates the 
basic physics of the interactions which would be present between a high 
voltage solar cell array and the plasma environment. Yet, unlike the solar 
array, the geometry is simple enough to be handled by a particle-pushing 

code. Although a ratio of ion to electron mass of 1:1 is commonly used, mass 
ratios of 100:1 or higher can be handled. This allows the simulation of 
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negative voltages present upon the central conductor. Thus, information 
gained- from an analysis of the disk and dielectric configuration- will be 
directly applicable to the design of large-scale., high-voltage solar. ceU 
arrays meant to be placed in low earth orbit. 


REVIEW OF THE-SIMULATION MODEL 


The plasma simulation code is a- further development of the 2-1/2 
dimensional program described previously (refs. 3 to 6). The cylindrical 
simulation region is divided into a numerical grid which is used to calculate 
the potential? and fields of the problem. The calculational grid is shown in 
figure 1. Each cell represents a ring of space due to rotational symmetry 
about- the 2-axis (boundary region 1) ... Boundary regions 2 and 3 represent the 
outer boundaries of the calculation. The outer boundaries are assumed to be 
removed far enough from the conducting disk so that potentials imposed upon 
the disk do not affect the boundaries. Thus, the plasma is maxwellian at the 
outer boundaries and the potential is set to zero. Particles are added to 
the simulation- along boundaries 3 and 4 according to the value of— the random 
thermal current. Particles are lost to the simulation whenever their orbits 
cross one of these two boundaries. Boundary region 4 represents the ground 
plate. The potential on boundary region 4 is set to zero and all particles 
intersecting this boundary are absorbed. The surface of the dielectric 
covering part of the conductor is represented by boundary region 5 (may have 
zero length). All particles striking boundary region 5 are absorbed. 
Electrons striking the dielectric may emit- secondaries with their number 
distribution given by fits to experimental data, presented by Haffner 

7 ) • The potential on the surface of the dielectric is determined by 
the following equation: 


®U) = ~ 4ira + -|«(0) + £$(az) (1) 


where 6 is the thickness of the dielectric, e is the dielectric constant, o 
is the charge per unit area, and 4>(0) is the potential of the conducting 
disk. This is just an infinite capacitor approximation which may not be 
valid near the edges of the dielectric surface. Boundary region 6 represents 
the exposed portion of the conducting disk which can be set at a desired 
potential with respect to the potential of the plasma exterior to the 
simulation region (zero) . Particles which enter boundary region 6 are 
absorbed and form the current drawn by the conductor from the-glasma. 

The simulation code is based upon standard PIC simulation techniques 
whenever possible. An overview of the program is shown in figure 2. The 
initial particles are randomly added to the simulated space with velocities 
chosen from a maxwellian distribution. This represents a nonequilibrium 
situation unless a potential of zero is specified on the conductor. The 
equations of motion are integrated using a second-order leap-frog method. To 
avoid the singularity present in the cylindrical equations of motion as r 
approaches zero, cartesian coordinates are used to compute particle 
movement. The coordinates are converted back into the cylindrical form when 


288 


ths move is completed. The particle mover advances particles through the 
simulation, in increments of the time step, value. The time step value is- set 
such that t|e fastest particles move less than 1 gjrid cell per time step 
(about 10 sec). Particles are added ... at the outer boundaries- from the 
surrounding undisturbed- plasma (assumed to have a maxwellian velocity 
distribution, usually of 10,000 deg K, about l_eV.). All- par tides present 
which have entered. .any of the sink regions are discarded. The remaining 
particle charges are. spread, over the grid using volume. weighting over the 
nearest grid cells. Successive over-relaxation with odd-even stepping is 
used to solve Poisson's equation - yielding a- self-consistent calculation of 
the electrostatic potential. The electric field can then be calculated by 
differentiating the potential. 

The total number of time steps simulated may range from about 4000 time 
steps for the plain disk configuration to 20,000 time steps for the disk and 
dielectric configuration. The average age of particles within the simulations 
for the 20x20 grid is about 100 time steps. This implies a complete cycling 
of particles in the simulation, space occurring about 200 times for a 20,000 
time step simulation. Hence, plasma parameters sucb as temperature and 
density are determined by the additions of particles on . the outside 
boundaries, not by the initial particle distribution. It is important to 
ensure that the program can maintain a constant density and temperature for a 
free plasma. aver many time steps. Simulation of a free plasma is accomplished 
by setting the radius of the conductor to zero (extending boundary region 4 
over the entire lower boundary) and reflecting any particles which cross 
boundary region 4. The potential on boundary region 4 is set to zero. A plot 
of the total number of particles present in the simulation region, versus time 
step number is shown in. figure 3. The initial number of particles (2000 per 
species) remains constant within statistical fluctuations for the entire 8000 
time steps computed using a 20x20 calculational grid size. The average 
kinetic energy per particle remains constant as shown in figure 4. Thus the 
simulation code has been shown to be able to simulate a free plasma for a 
large number of time steps without undergoing any nonphysical instabilities. 

This simulation program is applicable when the plasma Debye length is of 
about the same order of magnitude as the disk radius. Large plasma densities 
would cause large fluctuations in the collected current* since this 
calculation makes use of a small number of particles. The PIC assumptions 
made also imply that the Debye length is larger than one grid cell. With 
these restrictions and probe sizes ranging from 1.75 cm to 10.0 cm in radius 
and with electron and ioh temperatures of about 1 eV, plasma densities are 
then about 10,000 particles per cubic centimeter. This corresponds to 
environments encountered in the upper ionosphere or low earth orbit The 
Debye length under these conditions is about 5 cm. 


PLAIN DISK CALCULATIONS 


Approach to Equilibrium 

The plain disk or electrostatic probe configuration is the simplest 
configuration simulated. The rapid approach to equilibrium is demonstrated in 
figure 5. Voltage versus current density curves drawn by the disk are shown 
for simulation runs of various durations. The value for the current can be 
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the fact that the increased mobility of the electrons should charge the 
surface of the dielectric negatively. In order to be able to realistically 
simulate dielectric- charging effects and negative voltages for the disk and 
dielectric configuration, it is necessary to relax the PIC requirement of 1:1 
mass ratios. 

The disk only configuration provides an excellent test of the methods 
developed to handle larger mass ratios. This configuration is simple and 
gives simulation results with a minimum of CPU time. The expected collected 
current for each mass ratio can be found by simply scaling the^ posivive 
voltage result to the negative voltage simulation value. The diagnostics 
furnished are excellent. 

The usual way to handle larger mass ratios is to simply give the ions a 
larger mass value. This method has several disadvantages. Since the 
simulation procedure is not affected, it takes the same- CPU time per time 
step to simulate larger mass ratios as it would to simulate a 1:1 mass 
ratio. The distance each ion would move per time step decreases with the 
increasing value of the mass ratio. Thus, for reasonable mass ratios 
(hydrogen plasma, 1836:1), the simulation would have to be run- many times 
longer (at the same relative speed) to come to equilibrium.- Also the ion 
motion per time step might become so small as to be dominated by roundoff 
error. This generally limits PIC simulations to various small values of the 
mass ratio.- From these results, speculations are made as to the effect of 
realistic mass ratios. 

Another method to simulate large mass ratios can be developed. The 
relatively slow movement of the ions per time step indicates that it is a 
wasted effort to move them every time step. Instead, ions can be moved-once 
every n time steps. In this case n is chosen to be large enough so the ions 
move about the same average distances every n time steps as they would in a 
Simulation run with_a mass ratio of 1:1 for each time step. Thus ion and 
electron particle movement are restored to a relative parity at the expense 
of the additional assumption that the electric fields remain constant (at 
least in the average) for n time steps. The constant electric field 
assumption is clearly invalid at the beginning of a run. However, once a 
simulation reaches an equilibrium state, the assumption of a constant 
electric field should be justified. Once again it is postulated that if the 
simulation is run until an equilibrium state is reached, details about the 
approach to the equilibrium are lost - implying the equilibrium state is 
unique . 

The remaining problem is to determine n for a given mass ratio in such^a 
manner that the least amount of extra computing is required. The particle’s 
energy should remain constant as the mass ratio is varied. Therefore, the ion 
velocity is proportional to the square root of the inverse of its mass 
ratio. The time period over which the ions are moved then is just 
proportional to the square root of their mass ratio. The simulated mass ratio 
should be chosen such, that its square root is an integer. 

The above algorithm can be implemented with the PIC code by the 
following simple procedure. The form of the equations of motion describing 
the R coordinate which are solved by the particle mover are: 


291 


up 




m(V, 


r,Ql(P ~ At .| 

V \ 

*-S * 

^mr r l 

(2) 

old + At V 

,new 

(3) 


Since the energies of particles with differing mass ratios are constant at 
each point along the trajectory, the angular momentum term, l Z /mr , remains 
constant under variations of the ion mass ratio. The combined effects of 
substituting m'=n m and At ' =nAt just cancel in the equations of motion. A 
similar behavior is found for the equations of motion describing the 2 
coordinate. Thus, mass ratios of n :1 can be simulated simply by moving the 
ions only once every n time steps. 

Results obtained by using the direct :-d indirect methods to simulate 
large mass ratios have no effect on the current collected by the disk for 
5?®i tive voltages . The collected current is dominated by electrons: any 
difference in the ion mass is not seen. The best test of the larger mass 
ratio simulation methods is the simulation of negative voltages on the 
conductor. In this case the ions dominate the current and the electrons are 
excluded from the conductor. Results of simulations using several different 
mass ratios by both methods are shown in figure 8. The expected current which 
would be collected by the conductor decreases with increasing mass ratios. 
The general shape of the curxent-voltage curve remains intact for both 
approximations down to mass ratios of 4:1. For mass ratios of 100:1-, the 
straightforward method of computing mass ratio effects using a larger mass 
value begins to collect more current than the simulations which move the 
ions once per every 10 time steps. When a mass ratio of 1849:1 is attempted, 
by moving the ions once per 43 time steps, the current-voltage curve begins 
to bend over. This could be caused by a lack of statistics (30,000 time steps 
result in only 700 ion movements) or possibly the assumption of constant 
electric fields over the 100 time step interval is beginning to break down. 




THE DISK-AND DIELECTRIC CONFIGURATION 


aa. 
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The Approach to Equilibrium 

equilibrium situation of the disk and dielectric is much 
?h rd a r T t0 . ^ dentif y than that of the plain disk. The capacitance inherent in 
the dielectric causes the response time for charging to be large compared 
with the piasma frequency. Typically the dielectric may slowly charge or 
discharge during the simulation run, imitating an equilibrium situation in 

detail pH a ?iH£ S !!f fi the -i C Uf re S ValueS obtained * lt is important to take a 
detailed look at the available history information to determine whether a 

run reached an equilibrium state. For quick inspections, the value of the 
total charge collected by the dielectric is usually the most sensitive 
indicator of equilibrium. In practice it has been difficult to run 
current-voltage points long enough to arrive at an equilibrium value. 
Techniques have been developed to make the calculation come into an 
equilibrium state in a reasonable nu» *»r of time steps. 







, ma „ 0 , Way t0 re ? cb equilibrium is simply to run- the simulation long enough 
(many plasma- periods) until the equilibrium-state is fount!, This method was 
tried^ forthe 10x10 grid with SAG particles per species. A time history plot 
SJ*. 1 ,t ^ arg ? on the dielectric obtained from a 20,000 time step rim 

made at 10 volts is shown m figure 9. Beginning with an uncharged 

n.»n 1 t« t fh C ’ ^ arg ? a <- c,wm lat.es rapidly Cor about the first 15,000 time steps. 
l)ue to the uniformly decreasing value of the slope of the curve, it is 
difficult to determine the- onset of the equilibrium situation. 

h*,c«a T »f fi f? L !“ etlwd implemented to speed the approach to equilibrium is 
based upon the observation that the dielectric generally charges until its 
potential decreases to zero (the- voltage is lower than the first crossover 

potential. s ^ onda 7 em J ssl0n ^ >• If most of the charge needed to reach zero 
l aM f . duri " 8 the initialization, less time will be spent 
Jiff!™? 8 charg f dur J°* the ru n- Effects of loading the dielectric with 
different amounts of charge are also shown in figure 9. All runs come to the 

of me <w Ulll ^ riUm M whlch . su 88'* sts that the equilibrium is unique. This method 
of decreasing the time co reach equilibrium proves to be ineffective at 
higher voltages The dielectric must collect more particles to charge to an 

steps 'o^the ’cllLlation." 6 ^ **“ b0undarles duri “8 thf 2 M00 hm 

The relationship which causes the simulation to approach equilibrium 
slowly is contained in the boundary condition imposed upon the dielectric 
(equation 1). The leading factor of the thickness of the dielectric over the 

valie C of 1 i/?ft° nStant ^ fr USUally Sma11 (in the P rev ious runs we have used a 
valve of 1/28 as opposed to the experimental value- of 1/280 in order to keen 

mpanc *^1° J arge enough to be somewhat manageable). The small value of 6/e 
means that a large amount of simulation time will be used to accumulate this 
harge. The second method implemented to speed up the approach to equilibrium 

IfLcts of ntr0dUCt - 0n A/ artifiCiaUy high Values for fi/£ * The 4 physical 
effects of increasing 6/e can be thought of in one of several ways: 

1. a decrease in the capacitance associated with the dielectric 

2. an introduction of an artificial dielectric thickness 

Tho 3 ’ a w* ia }J y i®"**® 1 * 1 * the charge collected by the dielectric 

The result of a- 10 volt run with, a 10x10 grid of 500 particles per species is 

shown in figure 10 The time history graph of the total charge Collected by 

for a «lue ^ C (6/r)* S se^ ^ uilM,riU " is reached 500 time steps 

The artificial value of (6/e) must be chosen with care. If the value is 
the nnt*ent * equilibrium still will not be reached. If the value is too large, 
on h v P ° :“ J al ° n i. the dl ® lectric dramatically with the collection 8 0 f 

proportional to the thermal current leads to the following expression: 8 




(at)(n)(T) 1/z 


(4) 


where V is the voltage on the conductor, At is the time 
electron density, T is the electron temperature, and C 


step value, rj is the 
is the proportionality 
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constant. Substituting values of the other variables for 
run determines C (.4824). 


the above 10 volt 


The artificial value of (6/f,)‘‘ computed with the above equation assumes 
the number of particles present is large enough so that non physical effect- 
associated with- the granularity of chape collected on the dielcuic do 
occur. In practice, the value of (6/e) must, be kept small enough- that the 
change in the potential on the dielectric Surface per particle impact is onlv 

limi^i/l? 0110 ? ° f applied u ' tho conductor (1%) .^Thc ‘ uppej 

limit- of the value of (6/e) is then given as follows: 1 


(i\ * a 0.007 S(dr)V 

\c) mm —q— ( 5 ) 


where dr is the length of a grid cell and q is the charge per macroparticle. 


Comparison with Experiment 

f eSUltin 8 . current "VolUge curves for the disk and dielectric 

diele?tS l °radius Sh er ft%J igUre 1 V° r 3 disk radius of 1>75 cm and * 
r«nf^ ri Ji ad s of 8,75 cm * The simulation with a mass ratio of 1*1 
collects the same amount of current as the plain disk (the size of the 

exposed conductor is the same) for low- voltages. As the voltage on the 
disk UC surriunded Sed! th iS reached where th « current collected by the 

with positive^otentials near the^ondJc^ 
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Variations of the Secondary Electron Emission Yield Parameter 


scattered from the dielectric S surface m ner tile maxi,nu, “. number of electrons 
current voltage curves is shown in figure ™f 5332 iFtF'cSS *5 
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the increase in collected Current with reference to the plain- disk decreases 
with- increasing o . The increase in a lowers the energy of the fir>.t. 
crossover for secondary emission. Also, the number of secondaries release 
for energies larger than the first crossover point increases. Both effects 
create more secondary emission at lower voltages than would _ 

f und The sensitivity of the calculation to the secondary electron 

production might be sufficient to allow the measurement of the yield curve 
parameters for low voltages by fitting the response to experimental data. 

A detailed examination of the potential present along the top of the 
dielectric during a simulation run which collected more current whan that of 
the plain disk reveals what is happening. The initial condition shown in 
figure 13 for time step zero is a fully charged dielectric; the potential on 
the dielectric surface is zero. This would be realistic if the voltage o 
the conductor were increased slowly from a value low enough that the 
dielectric surface became fully charged. As the simulation proceeds, the 
potential upon the dielectric remains unstable, but positive for the tir 
6 000 time steps. After the equilibrium situation is reached, the potential 
remains constant to the termination of the run after 20,000 time steps. The 
equilibrium potential over the dielectric is uniformly positive and 
decreasing in value from the voltage on the exposed conductor to near zero at 
the outside edge of the dielectric. The potentials behave as if theyhav^ 
"snapped over" from their normal near zero level. This resembles the 
snapover phenomenon observed in experiments with solar cell arrays (ref. 2). 
The potential upon the cover slips for the array increased from values near 
zero below 100 volts to values about 5c volts less than that of the 
in ter con n ects for voltages greater than about 200 volts^ 


CONCLUDING REMARKS 


A cylindrical particle-in-cell plasma simulation code applicable to 
plasma densities encountered in low earth orbit is nearing completion. 
Results of the calculation of plasma coupling current for the plain disK are 
in agreement with experiment for positive voltages. Any deficiencies in e 
simulation of the disk and dielectric configuration are due solely to the 
interactions of the dielectric. Techniques have been developed which allow 
dielectric charging to occur at PIC time scales. The current-voltag 
characteristic curves are in qualitative agreement with experiment, 
indicating that for these voltage ranges charge sticking and secondary 
emission probably adequately describe the dielectric interactions with the 
ambient plasma. The amount of secondary emission from the dielectric at low 
voltages during the simulations needs to be reduced to match the experimental 
curreiit-voltage characteristic curve. Calculations using large ion to 
electron mass ratios are made possible by restricting both the grid size and 
particle number and simulating the mass ratio effects by moving the ions once 
for every n time steps. The square root of the desired ion to electron mass 
ratio determines the value of n. Further computational effort is required to 
extend the range of the disk and dielectric simulations to higher voltages 
and larger mass ratios. 
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Figure 2. - Overview of simulation code. 
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Figure 3. - Simulation of a free plasma - number of particles. 
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Figure 4. - Simulation of a free plasma - average kinetic energy 
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Figure 7. - Comparison with experiment - negative voltage 
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Figure 8.- Large mass ratio calculations. 
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Figure 10. - Approach to equilibrium - effects of i 
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- Disk and dielectric results. 
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Figure 12. ^ Variations in secondary emission yields. 
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Figure 13. - Snapover effect (potential snapshots) 



